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ABSTRACT 

A two-component hydrodynamic model is constructed of the global superfluid flow 
induced by two-component Ekman pumping during the recovery stage of a glitch. The 
model successfully accounts for the quasi-exponential recovery observed in pulsars like 
Vela and the "overshoot" observed in pulsars like the Crab. By fitting the model to 
high-resolution timing data, three important constitutive coefficients in bulk nuclear 
matter can be extracted: the shear viscosity, the mutual friction parameter, and the 
charged fluid fraction. The fitted coefficients for the Crab and Vela are compared with 
theoretical predictions for several equations of state, including the color-flavor locked 
and two-flavor color superconductor phases of quark matter. 
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1 INTRODUCTION 

Rotational glitches in radio pulsars unfold in two distinct 
stages: an impulsive spin-up event, followed by a quasi- 
exponential recovery towards steady spin down. The time- 
scales for these two stages are very different, suggesting that 
they involve different physics. The spin up is unresolved by 
existing radio timing observations but it is known to be less 
than ~ 40 s in obj ects like Vela, which h ave been mon- 
itored continuously (|McCulIoch et al.l 19871). The recovery 
typically lasts for days to weeks ( Wong et al.ll200ll ). 

The spin-up stage is generally attributed to a trans- 
fer of angular momentum from the neutron superfluid 
core to the crustal lattice vi a superfluid vortex unpin- 
ning IjAndreev fc Bashkinlll975r ). Recent data, which reveal 
that glitch sizes and waiting times follow power-law and 
Poissonian distribut ions respectively in individual pulsars 
(jMelatos et al.l2008l h suggest that the glitch trigger is an un- 
pinning avalanche in the superfluid vortex array, occurring 
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via a self-organized critical process 
200£ ) or a coherent noise process 
2009 ). To understand the collective nature of the trig- 
ger, one must synthesize a wealth of condensed matter 
physics at the microscopic level, including vortex instabil- 
ities, the distribution of pinning potentials, and the cou- 
pling of the superfluid to the charged fluid by entrainment. 
Separately, the spin-up event has been modeled hydrody- 
namically by averaging over the d i screte vortex micro physics 
(JGIampedakis fc Anderssonll2009l ; lsi"derv et al.ll2009l h In the 
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latter class of models, angular momentum is transferred 
from the superfluid to the charged fluid and hence the crust 
by mutual friction. 

The recovery stage is thought to reflect the restoration 
of sup erfluid-lattice corotation by viscous and/or magneti c 
forces (|Bavm et al]|l969l ; lBovnton et al.ll 19721 : lLohsenlll975T ). 
However, viscosity estimates based on electron-electron scat- 
tering and estimates of the magnet ic tension pre dict shorter 
recovery time-scales than observed (Easson 1979). Moreover, 
although mutual friction estimates based on electron scat- 
tering off neutron magnetic mome nts in vortex cor es are 
consistent with long recovery times (|Bavm et aI.lll969T ) scat- 
tering off vortex cores magnetized b y neutron-proton en- 
trainment (jAndreev fc Bashkinl 1 197a ) couples the proton- 
electron plasma to the neut ron superfluid on a time-scale 
much shorter than observed (jAlpar et alj|l984f ). In the lat- 
ter scenario, most of the fluid interior locks to the crust, 
and th e recovery stage is attributed to vortex repin ning and 
creep JAlpar et all 1 19931 . I1996J : ISiderv et ai]|2009l h All the 
above factors feed into the global hydrodynamics in ways 
that remain unclear. 

Three aspects of the glitch recovery process are truly 
puzzling from a physical perspective. First, if the r ecovery 
occurs hydrodynamically via the Ekman process (|Eassonl 
ll979l : lAbnev fc Epst ein 1996), one would naively expect the 
recovery time-scale to be the same for all glitches in a given 
pulsar, because the glitch amplitude and post-glitch per- 
turbation of the flow are small and hence the process is 
linear. The time-scale for linear Ekman pumping is inde- 
pendent of glitch amplitude and depends instead on con- 
stitutive properties such as viscosity and mutual friction, 
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which do not vary significantly during; the interval b etween 
glitches (JReiseneggerl Il993l ; lAbnev fc Epsteinl Il996h . Con- 
trary to expectations, however, the recovery time-scales ob- 
served in the Crab and Vela pulsars cover wide ranges 
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uOOll ; iDodson et al.ll2002r ) , a nd the same s eems to be true in 
less heavily studied pulsars (|Peraltall2007l h 

Second, glitch recovery cannot be parameterized by a 
single exponential decay. In almost all cases, at l east two ex- 
ponentials are required (|McCulloch et al.lll987), and up to 
four have been fitted to high resolution data ( Dodson et al.l 
I2002T ). The corresponding time-scales typically range from 
0.3 d to 300 d, suggesting that there are multiple physi- 
cal processes involved. This behavior is inconsistent with 
standard Ekman pumping: for example, the spin up of 
a stra tified viscous fluid in a cylinder involves one time- 
scale (jAbnev fc Epsteinl Il99q ). and the spin up of an un- 
stratified superflu id between parallel plates involves two 
(|Reiseneggerl Il993l ). The multiple time-scales are often at- 
tributed to the variation of vortex pinning strength (and 
hence vortex creep rate) with depth and have been modeled 
by dividi ng the star's moment of inertia into multiple com- 
pone nts (jAlpar et all 1 19931 Il996l ; ISedrakian fc Hairapetianl 
120021 ). 

Third, although many pulsars recover monotonically 
and quasi-exponentially like Vela, some pulsars do not. The 
Crab pulsar consistently "overshoots" during its recovery, 
decelerating below its stead y-state angular vel ocity before 
rising again asymptotically (|Wong et alj 120011 1. This sug- 
gests that the vortex unpinning event fails to redistribute 
angular momentum evenly throughout the star; differential 
rotation must persist between one or more internal compo- 
nents, which do not achieve corotation simultaneously dur- 
ing the recovery stage. 

In this paper, we consider an idealized hydrodynamic 
model for the recovery stage of pulsar glitches. Our model 
consists of a rigid outer crust containing a two-component 
superfluid. The viscous component (proton-electron plasma) 
spins up via Ekman pumping, whereby viscous stresses 
transfer angular momentum to the fluid in a boundary layer, 
which is then convected throughout the star by the Coriolis 



force (Greenspan & Howard 1963: !Reiseneggerlll993l : lEassonl 
ll979l : lAbnev fc Epsteinl 1996h . The inviscid component (neu- 
tron condensate) interacts with the viscous fluid compo- 
nent via the mutual friction force, which a rises from elec - 
tron scattering off magnetized vortex cores (|Mendelllll991f ). 
The important effects of vortex tension and pinning and 
macroscopic entrainment are neglected to keep the prob- 
lem tractable analytically. In contrast to previous stud- 
ies, the crust has a finite mass and responds to the vis- 
cous torque applied by the interior fluid. This back-reaction 
in turn modifies the flow field and we calculate this self- 
consistently. The magnetic dipole torque is neglected, be- 
cause the recovery stage is much shorter than the electro- 
magnetic spin-down time. Our approach differs from pre- 
vious studies of the recovery stage, which incorporate vor- 
tex pinning and creep but do not solve self-consistently for 
the global two-fluid flow pattern produced by Ekman pump- 
ing ISedrakian fc HairapetiadlJffloll ; lAlpar et al]|l984ll993l . 
ll996T ). The global flow pattern and viscous back-reaction are 
built into a hydrodynamic model self-consistently for the 
first time, yielding a tool that can be used to investigate a 



variety of pictures for pulsar interiors, as we discuss below. 
We do not assume a priori knowledge of any constitutive 
coefficients, leaving them completely free to be determined 
by observations. Thus, we seek to determine what aspects 
of the post-glitch relaxation a self-consistent hydrodynamic 
model can and cannot explain, in order to clarify what ele- 
ments of non-hydrodynamic physics are absolutely required 
by the data. 

In 321 we present an analytic model to describe the ro- 
tational evolution of the crust during the recovery stage. 
The model is based o n a recent analytic solution by 
IVan Eysden fc Melatosl (|2010f ) for the spin-up problem of 
a two-component superfluid in spherical Couette geometry. 
In this paper, we specialize to the case where there is no in- 
ner core, saving the more general problem for a future study. 
The output of the model can be compared directly against 
high-resolution radio timing data. It is fitted to the quasi- 
exponential recovery of Vela glitches in §3] and to the "over- 
shoot" recovery in Crab glitches in Sj3] In both cases, the rel- 
evant superfluidity coefficients are extracted. Finally, in <}5] 
we compare the fitted values of the superfluidity coefficients 
with the oretical predictions for a s t andard 1 5'o ne utron su- 
perfluid (|Cutler fc Lindbloml 1 19871 : iMendelll 1 199 if ) and for 
various dissipation channels in an exotic str ange-quark su- 
perfluid, e.g. in a color flavor locked phase (|Madsenl |2000| ; 
iMannarelli et al.l 120085 : I Afford et ai]|2009l '). There is cur- 
rently a flowering of interest in dissipative transport pro- 
cesses in bulk nuclear matter, as the foregoing references in- 
dicate. Quantitative glitch recovery studies, especially when 
performed on pulsars with different ages, offer a promis- 
ing way to measure nuclear transport coefficients experien- 
tially in the many-body, MeV regime, which is inaccessible 
at present in terrestrial laboratories. 



2 SPIN UP OF A SPHERICAL STAR: 
TWO-FLUID THEORY 

In the absence of a satisfactory microscopic explanation of 
the glitch trigger, we consider the hydrodynamic response 
of the stellar interior to the following set of idealized ini- 
tial conditions. Consider a thin, spherical shell of radius R, 
representing the solid crust of the star, which rotates at an- 
gular velocity Q just before the glitch (t = 0~). Suppose 
the shell contains a two-component superfluid, whose invis- 
cid (Bose-Einstein-condensed neutrons) and viscous (uncon- 
densed neutrons, protons and other charged species) compo- 
nents rotate rigidly but differentially with angular velocities 
f2 s o and O„o 7^ fl s o respectively at t = _ . Immediately af- 
ter the glitch, the crust accelerates instantaneously to reach 
an angular velocity f2 + AQ at t — + . Subsequently, at 
t > 0, the angular velocity of the crust evolves according to 
O + AQf(t) [with /(0 + ) = l], while the two superfluid com- 
ponents participate in coupled Ekman pumping, with the 
coupling provided by mutual friction. During the Ekman 
process, neither component rotates rigidly, unl ike in other 
"body-averaged" treatments (jSiderv et al.ll2009r i. At its con- 
clusion, however, the components tend toward corotation. 

The above initial conditions, although idealized, are 
consistent with the spirit of the superfluid vortex unpin- 
ning paradigm, in the sense that the inviscid component 
leads the other components before the glitch and transfers 
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part of its excess angular momentum almost instantaneously 
(e.g. via Kelvin wa ves) to the crust during the vortex un - 
pinning avalanche (JAlpar et al J 1 1984 iMelatos et al]|2008h . 
The absence of a "reservoir effect" (glitch size ex waiting 
time since previous glitch) in the data guarantees that the 
avalanche event nullifies only a small fraction of the accu- 
mulated differential rotation a nd leaves £l s o ^ il n o at t — + 
(jMelatos fc Warszaw ski 2009). (The instantaneous decrease 
in Sl s at t — accompanying the avalanche can be absorbed 
into fi s o without loss of generality.) By the same token, the 
final conditions (corotation between all components) are less 
realistic; they fail to acknowledge that differential rotation 
must persist at all times to avoid the reservoir effect. This 
flaw is shared by all published hydrodynamic models, which 
elect not to describe the stochastic vortex repinning process. 
In what follows, we adopt cylindrical coordinates 
(r, <f>, z) and work in the non-inertial frame rotating with the 
pre-glitch angular velocity Qk of the crust, where k denotes 
the unit vector along the z-axis. 



2.1 Interior flow 

The flow in the interior of the pulsar is governed by the two- 
fluid Hall-Vinen-Bekharevich-Khalatnikov (HVBK) equa- 
tions. These equations can be linearised on the basis that 
the glitch amplitude is small (Afi -C O). In the rotating 
frame, the incompressib le, linearized HVBK equations take 
the dimensionless form (jPeralta et al.l l2008) 

E i/2 ^_ + 2kxVn = _ Vpn + EV 2 v n + p s F, (1) 

OT 



E 1 ' 2 -^- + 2k x v s = - Vp s - p n F , 
or 

V • v„ = , 

V • v 3 = . 



(2) 
(3) 

(4) 



The components ar e coupled by the mutual friction force 
JHall fc Vinenlll956l ) 

F = Bk x [fc x (v n - Vs)} + B'k x (v n - v s ) , (5) 

where we assume the flow is laminar [cf. IMelatos fc Peraltal 
(|2007l )] and F takes the anisotropic Hall-Vinen form rather 
than the isotr opic Gorter-Me l link f o rm [i .e. , there is no vor- 
tex ta ngle; cf. iPeralta et all (|2005l . |2006| ); lAndersson et al.1 
(J2007T) ] . The symbols «„ )S and p„, s denote the velocity and 
pressure of the inviscid (s) and viscous (n) components. B 
and B' are dimensionless parameters which we hope to mea- 
sure by applying the spin-up model to pulsar timing data. 
The velocity and pressure scales are chosen to be RAfl and 
pQR 2 AQ respectively. The mass densities of the viscous and 
inviscid components, p n and p s , are scaled to the total den- 
sity, p, so that we have 



p n + ps = 1 ■ 



(6) 



To scale the time coordina te, we define the Ekman time 
(jGreenspan fc Howa rd 1963) 



t = E 1/2 nt 

in terms of the Ekman number 

/' 



E 



PnR 2 ^ ' 



(7) 

(8) 



where p is the shear viscosity. In a neutron star, E is ver y 
small, with E < 10~ 10 typically JMelatos fc Peraltal [20071 ). 
Note that the centripetal force terms are absorbed into the 
definitions of the pressures, as the flow is incompressible. 

Equations ©-(SI neglect entrainment, whereby the 
flow of one fluid imparts momentum to the other and 
vice vers a via a quantum mechanical current-current in- 
teraction ( Andersson et al.ll2006l ; lAndersson fc Comerll2008l ; 
ISiderv et al.ll2009l ). Entrainment is expected to play an im- 
portant role in the dynamics of neutron star interiors, but 
we neglect it in this paper to keep an already difficult 
problem analytically tractable. [Superfluid spherical Cou- 
ette flow has never been solved analy tically before and w as 
treated numerically only recently; see lPeralta et al.l |2008l ).] 
Similarly, magnetic fields are also neglected, even though 
Ivan Hoven fc Levinl (|2008l ) showed that they are impor- 
tant in vortex dynamics, e.g. suppressing the Donnelly- 
Glaberson instability. Vortex tension is neglected on the ba- 
sis that the inter- vortex spacing is small, although the exact 
nature and str ength of vo rtex pinning to the crust is still 
being debated (|Linkll2009l ). 

The presence of the mutual friction force in CO)-© 
introduces new physics into the traditional Ekman pump- 
ing process. In addition to the classical Ekman time-scale 
E^ 1 ' 2 ^ 1 , the mutual friction introduces a second time- 
scale, B _1 S1 _1 , characterized by the coupling strength. The 
ordering of these time-scales governs the dynamics of the 
spin-up flow. For B ~ 1, the inviscid and viscous fluids are 
locked together; differential rotation is removed by mutual 
friction over a few rotation periods. For B ~ E 1 ' 2 , the spin- 
up time is a combination of E^ 1 ' 2 ^ 1 and B^ 1 ^ 1 ; the 
viscous fluid spins up via Ekman pumping, while "drag- 
ging" the inviscid component along via the mutual fric- 
tion. For B -C E 1 ' 2 , Ekman pumping proceeds for the 
viscous fluid uninhibited by mutual friction over the time- 
scale E" 1 ' 2 ^!" 1 , while the inviscid component is brought 
into corotation over the much longer time-scale B -1 f2 -1 . 

In helium II, the mutual friction arises when excited 
states scatter off the vortex lines, and B and B' are of or- 
der unity. In neutron stars, several mechanisms give rise to 
mutual friction, including: (1) electron scattering off the 
neutron magnetic moment s in vortex cores , with charac- 
teristic time-scale ~ lyr JBavm et al.l 1 19691 ); (2) electron 
scattering off entrained protons, which magnetize the vor- 
tex c ores, with time-scale ~ Is |Afpar et al.lll984l ; iMendelll 
Il99ll ); and (3) excited states scattering off vortices like in 
Helium II. Processes (1) and (2) bring the neutron superfluid 
into cor otation with the p lasma via electromagnetic inter- 
actions (|Alpar et al.lll984T ). Uncharged particle species (e.g. 
excited states of the neutron condensate or exotic particle 
species) couple via (3). Terrestrial experiments on liquid he- 
lium reveal that ~ 10% of the fluid is in excited states even 
at absolute zero, well above the ideal Bose gas fraction, be- 
cause of the non- ideal nature of the molecular forces. Similar 
non-ideal behavior for the strong nuclear force is not ruled 
out at present. 

Transport coefficients such as viscosity and mutual fric- 
tion depend on depth (via the density and temperature, 
see ij5]), and neutron stars are strongly stratified. Hence, 
by treating E, B and B' as uniform in this paper (to 
keep the problem tractable), we are obliged to interpret 
the values generated by fitting the model to data (see ij3] 
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and Q as body-averaged effective values. The average is 
taken over regions of the star with different compositions 
(e.g. exotic quark matter), densities, and temperatures. Im- 
portantly, the average depends sensitively on the depth 
to which Ekman pumping extends; it is reduced dramati- 
cally by compressibility and stratific ation (|Abnev fc Epstein! 
1 19961 ; IVan Eysden fc Melatosl I2010T ), crucial physics which 
we do not consider here. 

The general analytic solution to CH) — (SI) m 
spherical Couette geometry was derived recently by 
IVan Eysden fc Melatosl (120100. generalizing the boundary 
layer analysis of iGreenspan fc Howard! (|1963l ). In this paper 
we restrict our attention to the simple situation where there 
is no inner core and specialize to the regime B,B' -C 1 
(which we verify a posteriori in several examples in O 
and jj4l) pertainin g to glitching pulsars (|Mendel]| Il99ll ; 
ISiderv et al.l 120091 ). For the general solution with no inner 
cor e, the reader is present e d the Appendix. As discussed 
by IVan Eysden fc Melatosl (|2010T l. the solution to JTJ-© 
is expressed most neatly in terms of the total mass flux, 
defined as 

v - p„v„ + p s v s . (9) 

The azimuthal component of the total mass flux is given by 

ru)+ {r)ui- (r) 



v^{r,r) = 



dr'f(r') 



0[u,+(r)-u,-(r)]j 



(10) 



rQ. n0 Lj + (r)uj-(r) 



P[u+(r)-u-(r)] 

rQ 



[« 



_(r)r 



[w_(r)e 



-(Or 



■ Wt 



■u+(r)e' 



,_(r)r] 



w+(r) - w_(r) 

In (|10p . Qo and Q„q denote the initial angular velocities 
of the total mass flux and viscous component in the rotat- 
ing frame, scaled to Af2. The radius-dependent time-scales 
u)^ 1 are a mixture of the classical Ekman time-scale and the 
superfluid mutual friction coupling time-scale as discussed 
above, with 



u±(r) = -- 



± 



P+(l- r 



p+(l-r 



-3/4 



-3/4 



- Pp n (l - r 



-3/4 



(11) 
1/2 



and 

P = BE~ 1/2 . 



(12) 

Note that B' no longer appears in the equations in the 
regime B, B' <C 1. The boundary conditions leading to (I10[) 
are that the viscous component v n satisfies no penetration 
and no slip at Z — ±(1— r ) , while the inviscid component 
v s satisfies no penetration. 

We quote only the result for the interior azimuthal flow 
in ()10[) . The full solution also consists of the radial and ver- 
tical flows describ ing the Ekman flow, as well a s boundary 
layer corrections (|Van Eysden fc Melatosl 120101 ). However, 
the latter elements do not appear explicitly in the theory 
of glitch recovery. 



2.2 Back reaction torque on crust 

Equation (|10p gives the solution of (|T|)— (JJJ) in terms of a 
general boundary condition /(r). To solve the glitch recov- 
ery problem self-consistently, we must determine the time 
evolution of /(r) due to the hydrodynamic torque on the 
rigid crust. This is done by integrating the viscous stress 
te nsor over the stellar surface z = ±(1 — r 2 ) 1 ' 2 . The result 
is (|Van Eysden fc Melatodl201Ch 



df(r) 
dr 



15K 



drr (l 



2\i/2 dv^{r,r) 
" I dr 



(13) 



where K is the ratio of the moment of inertia of the to- 
tal fluid to the moment of inertia of the crust. Equations 
(|10|l and (|13|l constitute a closed pair of integro-differential 
equations for the unknowns v^ir^) and f(r). Substituting 
(|10|) into (113[) . we derive an integral equation for the spin 
evolution of the crust, viz. 



f(r) = -p n K / dr' [g A (r - r') + g B (r - r')] f(r') (14) 
Jo 

+ p n K [g A {r)n nil + g B (T)Q ] + 1 , 



where we define 
9 (r) = y 



dr 



b 



„(r)T 



.(r)-r 



li 



(l_r2) 1/4 [« + (r)-w_(r)] 



9 (j) 



P 



/ dr'p A (r') 
Jo 



(15) 



(16) 



It is straightforward to solve (|14[) for /(r) numerically, by 
guessing an initial trial function for /(r) (e.g., exponential), 
substituting it into the r ight-hand side o f l|14p . updating 
/(t) via underrelaxation (jPress et al.ll2002J '). and iterating. 

The parameter K can be interpreted in a number of 
ways, depending on the glitch microphysics. It can repre- 
sent the moment of inertia of the crust (i.e., ionic lattice 
alone), or the total moment of inertia of the crust and all 
components of the star coupled to it on a very short time- 
scale. Broadly speaking, if the crust is stro ngly coupled 
magnetically to th e proton-electron plasma (|Bavm et al.l 
1 19691 ; lEassonlll979l ) and subsequently to other components 
of the st ar, via electron sca ttering off magnetized neutron 
vortices (lAlpar et al.l ll984T ) or vortex-fluxoid interactions 
l|van Hoven fc L evin 2008), one has K ~ 1. The viscous fluid 
component then represents excited neutron states (non-ideal 
internuclear forces) or other neutral exotic particle species 
decoupled magnetically from the crust (see i]4.5[) . Alterna- 
tively, if K refers just to the ionic lattice, one has K > 50 
(see p. El 

An advantage of the above model is that is predicts the 
behavior of an observable, namely Q + Af2/(r), the post- 
glitch angular velocity of the crust. The model has six free 
parameters: p n , K, B, E, fi„o and Qo- 



1 In a more general theory where the magnetic dipole spin-down 
torque is included, it would act on the crust and any components 
strongly coupled to it. 
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2.3 Fitting to radio timing data 

It is customary to fit the post-glitch crustal frequency v g (t) 
empirically using a sum of N exponentials, with amplitudes 
and e-folding times Av n and t n respectively, i.e., 



where /(oo) is defined in (|20|) and C is given by 



Vg(t) = VgO + Av p + \J Al 



-t/tr, 



(17) 



where Av v is the permanent part of the spin up, v g o = 
v g (0~) is the spin freq uency just before the gl i tch, and one 
has N ^ 4 typically (IMcCulloch et all I1987J ; lAlpar et all 



1 19931 . H996; W ong et al 



20011 s ). Normalizing (fTT)) to the ob- 



served frequency jump 



Av = AQ./2-k 



N 
Av p + 2_^ Al/n 



(18) 



immediately after the glitch, and writing t in terms of the 
Ekman time ([8]), we obtain 

A N A 

Us(r) = ^+J2=^ exp{-T/lE 1/2 2w g0 t n ]} . (19) 

n = l 

Our goal is to adjust the six model parameters p n , K, B, 
E, Q n o and flo until the analytic solution f(r) from (|14[) 
matches the observed behavior fobsij) as accurately as pos- 
sible. 

Post-glitch timing typically yields ~ 10 2 data points (~ 
1 per day) during the recovery stage, with formal residuals 
< 1%. Hence the six model parameters are overconstrained 
in principle. In practice, the comparison of fobs(r) and f(r) 
is done by eye, and the idealized theory in t]2.1l and ^2.21 is 
best matched to the approximate observed behavior given 
by dT7|. 

Two conditions on the model parameters can be read 
off the data straight away. First, con servation of angular mo- 
ment um in the steady state implies (|Van Evsden fc Melatosl 
l2010h 

Av p 1 + KQo 

a^ = /(oo) = TTF- 



(20) 



Second, the initial slope 

ijVan Evsden fc Melatosll2010h 



/(0 + ) 



satisfies 



Ell A "n/t n 



£ x// /(o + ) 



20 
' 7 



p n KE 



1/2 



2-KVgoAv 

Equations (|20p and (|2ip provide two conditions on the model 
variables (on the right-hand sides) in terms of the mea- 
sured quantities Av p , Av„, t n , v g o (on the left-hand sides). 
In general, the remaining four conditions must be deter- 
mined through trial and error by matching fobsij) and f(r) 
by eye. 

2.4 An approximate solution involving dual 
exponentials 

If the superfluidity coefficients satisfy E 1/2 -C B, B' -C 1 
(which we verify a posteriori in several examples in Sj3]and 
and K ^> 1, then (|14|l has the approximate solution 

f(r) = [1 - C - /(oo)] exp[-(20/7K(l + K)r] (22) 

+ Cexp(-B£~ 1/2 T) + /(oo) , 



C = 



20 Pn K (Q - Q n0 ) 
7BE- 1 / 2 - 20p n K 



(23) 



This approximation matches the exact numerical solution to 
s£ 1% for K > 10 3 . Comparing (JT5J and p2)>. we can read 
off the model parameters directly in terms of t\ and t-^: 



B = (2-KVgoUy 



20 



— E 1/ 'p n (l + K) = (2nvg Q t 2 ) 



(24) 
(25) 



Equations Q24p and (|25p provide two more observational con- 
straints in addition to (|20[) to (121[) for the model parame- 
ters. The model is therefore underconstrained; we are free to 
choose two parameters. For any given K and p n , say, equa- 
tions IpO). (f2T|) . ([24)l and (J25J) determine B, E, fl n0 and O . 
Note that a second, equally valid set of model param- 
eters can be extracted by taking the time-scales the other 
way around, i.e. by swapping t\ and £2 in (|24|) and (|25[) . 
The approximate solution also provides an excellent initial 
trial function when iterating (|14[) . Of course, the full tim- 
ing ephemeris contains more information than equation (|17l) 
with N = 2,so the model is actually overconstrained in prin- 
ciple, as noted above. Nevertheless, the quality of the data 
and realism of the model are such that two of the param- 
eters are free in practice in many objects. If we insist that 
K and p n (say) do not change from one glitch to the next 
in an individual pulsar, then we can constrain all six model 
parameters uniquely. This is done in ij3]and S|4] 



3 VELA 

3.1 Data 

Since 1969, Vela has been seen to glitch a total of 17 
times, comprising 15 macro-glitches (Av ~ 10 -e V 9 o) and 
two micro-glitches (Av ~ 10 _8 i' g o). The modified Ju- 
lian date and timing parameters of each event are listed 
in Table [l] The first four glitches were discovered in 
data collected by the Deep Space Network (ICordes et al 



1988; [Radhakrishnan fc Mancheste r 1969; Manchester et al 



ttno) ■ (21) 



1976T I. Since then, Vela has been monitored almost con- 
tinuously by two radio telescopes: the 14-m antenna at 
the University of Tasmania's Mt Pleasant Observatory 
(5 hours p er day from 1981 to 1987, 1 8 hours per day 
thereafter) (IMcCulloch et al.l Il983l . 1 19871 : IMcCulloch! ll996J : 
iDodson et al.ll2002l , 120071 ), and the 26-m antenna at the Har- 
tebees thoek Radio Observatory (11.5 hours per day s ince 
1984) (|Flanaganlll9"90l . ll996l : lBuchner fc Flanaganll2008l ). In 
their present configurations, the telescopes record all four 
Stokes parameters at 635/950/1390 MHz and 1.644/2.326 
GHz respectively. To date they have each detected a to- 
tal of nine glitches. Timing p arameters have not yet bee n 
published for the latest event (jBuchner fc Flanagan! 120081 ") . 
which does not appear in Table [l] 

Continuous timing data is crucial for the work in this 
paper, because an accurate measurement of the glitch epoch 
feeds through into accurate measurements of the recovery 
time -scales t n , the shortest of which can be a fraction of a 
day (|Dodson et al.ll2002l , [20071 '). All glitches for which data 
were recorded during the event itself are marked with an 
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Table 1. Timing parameters for large glitches (Au ^ l^tHz) in the Vela pulsar (u g o = 11.2 Hz). An asterisk in the CM column denotes continuous monitc 



Glitch Date 



MJD CM 



*4 



t 3 



t-2 



Av± Al^3 



Ai^2 Av\ 
/iHz 



Ai/„ 



Av Ref. 



1 


28-Fcb-1969 


40280 


2 


29-Aug-1971 


41192 


3 


09-Sep-1975 


42664 


1 


13-Jul-1978 


43693 


5 


10-Oct-1981 


44888 




10-Oct-1981 


44889 


6 


10-Aug-1982 


45192 




10-Aug-1982 


45192 


7 


12-Jul-1985 
12-Jul-1985 


46258 


8 


24-Dec-1988 


47519 




24-Dec-1988 


47520 


9 


20-Jul-1991 
20-Jul-1991 


48458 


10 


26-Jul-1994 
26-Jul-1994 


49560 


11 


27-Aug-1994 
27-Aug-1994 


49592 


12 


13-Oct-1996 


50370 


13 


16-Jan-2000 


51559 


11 


07-Jul-2004 


53193 





10 


120 




1 


91 




1 


35 




6 


75 




6 


11 




1.6 


233 




3 


21.5 




3.2 


60 




6.5 


332 




6.8 




0.4 


4 


96 


0.73 


6.97 


707 


0.56 


5.94 


251 


0.59 


4.9 


49 



0.108 
0.092 
0.255 
0.317 



0.0008 
0.0007 



0.53 
0.23 



1.59 



3.29 
2.1 



15 

916 
19 

26.14 



0.02 

54 



0.31 
0.21 



0.052 
0.036 

0.01 
0.083 

0.01 
0.092 
0.057 

0.23 
0.066 
0.165 
0.086 
0.083 
0.169 
0.152 



0.021 



0.193 
0.13 



0.467 

0.30 

0.079 

0.389 

0.024 

2.26 

0.126 

0.79 

2.76 

0.376 
6.74 
2.84 

0.231 



25.6 
22.6 
22.2 
33.8 
12.7 
10.5 
22.8 
22.0 
15.1 

19.7 
13.3 
27.1 



0.027 2.1 
0.032 

14.8 9.1 

0.236 34.5 

0.16 22.8 



26.2 
22.9 
22.2 
34.3 
12.7 
12.8 
23.0 
23.1 
17.9 

20.2 
20.2 
30.3 



9.6 9.6 



2.2 

23.9 
35 

77.3 



[1] (Cordes et al. 1988) . [2] (McCulloch et al. 1987) , [3] (Flanagan 1996) . [4] (Flanagan 1990) . [5] (McCulloch 1996) 
[6] (Wang et al. 20001 . [7] (Dodson et al. 2002) . [8] (Dodson et al. 2007) 



asterisk in the fourth column of Table \T\ The timing pa- 
rameters in Table [1] are defined in terms of equation (|17[) . 
Note that the Hartebeesthoek group fits for the derivatives 
Az> n rather than Av n ; the former are converted to the lat- 
ter in Table [T] Some events were observed simultaneously 
at Tas mania and Hart ebeesthoek, e.g., the 1988 Christmas 
glitch (|Flanaganll 19901 ) . As the efforts to continuously mon- 
itor Vela have intensified, the timing parameters have been 
determined ever more accurately. A third, short time-scale 
£3 has been discernible since 1988, together with a fourth, 
even shorter time-scale £4 since 2000. The exceptions are 
the "double glitch" of 1994, where no discernible relaxation 
was observed in the month between the two events (hence 
empty entries in Table [1} , and the 1996 glitch, which was 
not observed at either Mt Pleasant or Hartebeesthoek. Ta- 
ble Q] also presents two different sets of timing parameters 
for the 1988 Christmas glitch. Both Mt Pleasant and Har- 
tebeesthoek were monitoring continuously at the time and 
fitted similar values of Az/2 and Az/3 but very different values 
of A1/1 and Av p . This discrepancy appear s to arise fr o m the 
tail of the rec overy, whi c h was tracked bv lMcCullochl (| 19961 ) 
for 707 d and lFlanaganl (|l996T > for 96 d. 

It is important to bear in mind that, just because the 
empirical timing model (|17[1 adequately fits the data, it 
does not mean that this functional form is representative 
of the underlying physics of the recovery stage. As more 
time-scales are revealed by higher resolution data, the pos- 
sibility grows that the true functional form for the recovery 
may be something other than a sum of exponentials, even 
if it is well approximated by the latter. One such possibil- 
ity is a "non-linear" decay term that ties the intermediate 
and lon g-term relaxation to weak and superweak pinning 
regimes (|Alpar et al.lll993l ). A second possibility is explored 
in the remainder of this section, where we show that (i) the 




Figure 1. The 1985 Vela glitch, as fitted by a one-component 
spin- up model ( H3.20 . showing the normalized spin frequency 
f(t) as a function o f time (in d). The data collected by 
iMcCulloch et al.l |l987J) are displayed as a heavy black curve. The 
lighter curves correspond to the theoretical model with K = 10~ 2 
(blue, top) and K = 10 2 (red, bottom). The dashed curve corre- 
sponds to a pure exponential spin down matched to the observed 
slope at t = 0; it does not coincide with either the data or the 
theory. The dotted curve is the steady-state spin frequency /(co). 



simplest hydrodynamic model predicts that v g (t) is not ex- 
ponential, and (ii) the complete, two-component superfluid 
model reproduces (|17|l in a natural way for N ^ 3. 



3.2 One-component viscous fluid 

We begin by investigating the classical case where there is no 
superfluid component. The explicit forms of equations (IIOII 
(|16|l in this limit are presented in SjB] The equations involve 
only one free parameter, K. The classical model is compared 
to observational data in Figure [TJ using the 1985 glitch as a 
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test case (|McCulloch et alJll987l ). The data themselves are 
plotted as a thick black curve, based on the timing parame- 
ters in Table [l] A pure exponential decay which matches the 
initial slope and steady state is plotted as a dashed curve. 
The remaining curves correspond to K — 1CP 2 (blue, top) 
and K = 10 2 (red, bottom). 

Figure Q] demonstrates that a spherical crust experienc- 
ing a viscous Ekman torque does not spin down exponen- 
tially, even if the spin- up is linear (Af2 <C fl); the red and 
blue curves and the dashed curve in Figure [l]do not coincide. 
This non-exponential behaviour occurs because the Ekman 
torque varies with latitude in a sphere. Hence, even in the 
simplest model, there is a natural hydrodynamic explana- 
tion for multiple time-scales, namely the Ekman time-scale 
associated with the "ring" of crust at each latitude. 

Figure \T\ demonstrates that there is no value of K which 
fits the observational data for the 1985 glitch in the simple 
model. The same holds true for all the other events in Ta- 
ble [T] Therefore, we must look to more detailed models to 
explain the form of v g (t) observed during the recovery stage. 



3.3 Two-component superfluid 

Next we apply our complete two-fluid model [equation (|14[) ] 
to the 1985 test case. In Figure^ w e present two dist i nct fit s 
to the timing solution measured bv lMcCulloch et al.l (|1987l ). 
The model parameters are given in the caption. The two fits 
correspond to (a) a heavy crust, with K — 1 and BE -1 ' 2 = 
23.6, and (b) a light crust, with K = 50 and BE~ 1/2 = 739 
with pn = 0.1 in both cases. In both fits, the theoretical 
curve is drawn in blue and the data is drawn in red. Also 
plotted is the approximate solution given by (|22[) . which is 
valid in the regime E 1/2 < B,B' < 1 and K » 1. We 
see that (|22|) approximates the exact solution more closely 
in Figure 2(b), where BE -1 ' 2 and K are larger, than in 
Figure 2(a). 

Interestingly, in the regime E ' <C B,B' <C 1 and 
K S> 1, /(t) reduces to the sum of two exponentials. Yet, 
for a single viscous fluid, /(r) is the convolution of many 
exponential time-scales at different latitudes (see H3.2p . Why, 
then, does the multi-exponential behavior go away, when the 
two-component superfluid still contains a substantial viscous 
component? The reason is that the mutual friction force and 
the relatively large fluid inertia dominate the flow dynamics 
in the above regime. As a result, the latitudinal variation of 
time-scales associated with the viscous torque is insignificant 
in relative terms. 

Figures [2ja) and [2jb) imply that there is a degeneracy 
in the fitting parameters for this glitch. For increasing values 
of K, a collection of accurate fits may be obtained. Indeed, 
it transpires that an infinite, one-parameter family of valid 
fits exists. This can easily be understood in terms of the 
approximate solution. As discussed in £12.41 when only two 
time-scales are resolved in the timing data, / t s (r) can be 
compared directly to (|22p to yield (|24|l and (|25|l . Combining 
(1241) and (1251). we arrive at the relation 



BE 



-1/2 = ( 20t 2 



V 7ti 



Pn(l + K). 



(26) 




p„(l+K) 

Figure 3. Superfluidity coefficients BE -1 ' 2 and Pn(l + K) in- 
ferred from the 1985 Vela glitch. The points are obtained by fitting 
the theoretical solution / (t) from l|14| l to the data by eye. The 
open boxes correspond to the fits in Figures [2f a) and (b). The 
filled circles indicate five additional fits. The two curves are plots 
of the relation (126 I t obtained by fitting two exponentials to the 
data, with (ti,t 2 ) = (6.5, 332) (top) and (332, 6.5) (bottom) (t\ 
and i 2 in units of d). 



1.000 
0.995 
0.990 



which predicts a direct proportionality between the su- 
perfluidity coefficients BE" 1 ' 2 and p n (l + K), with slope 
(20t 2 )/(7ti). 



20 40 60 80 100 120 

t (days) 

Figure 4. The 1988 Vela glitch, as fitted by a two-component 
spin-up model S13.3I showing the normalized spin frequency fit) 
as a function of time (in d). The blue curve represents /(£) for 
B = 5.0 x 10~ 9 E = 3.0 5 x lO " 15 , p n = 0.01, K = 53, tt = 
f2 ra o = 0.97. The lFlanaganl (1990) data / j, s (t) are graphed in red. 
The dotted curve is the steady-state spin frequency /(co). The 
agreement is good. 



In Figure [3] we plot (|26p as a straight line on the 
BE~ 1 ' 2 -p n (l + K) plane on a log-log scale. We also plot 
the parameters obtained by fitting /(r) to the 1985 data 
by eye, including the fits in Figures 2(a) and 2(b) (open 
boxes) and five additional fits (filled circles). The param- 
eters extracted from these fits are quoted in Appendix C. 
We find that (|26|) is consistent with the by-eye fits, even 
in the regime BE -1 ' 2 ~ 1 where the approximate solution 
nominally breaks down. 

The lower line plotted in Figure [3] is obtained by swap- 
ping t\ and £ 2 in (|24[) and ([25), as discussed in £12.41 By-eye 
fits are harder to achieve along this second line, because the 
iterative numerical solution to (|14jl is unstable. 

Does the two-component model perform equally well 
for glitches with three time-scales? In F igure [4] w e prese nt a 
fit to the timing solution measured by iFlanaganl (|19901 ) for 
the 1988 Christmas glitch. The model reproduces the data 
admirably. Moreover, the degeneracy expressed by (|26|) is no 
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400 600 

1 (days) 




1 (days) 



Figure 2. The 1985 Vela glitch, as fitted by a two-component spin-up model f i|3.3l) . showing the normalized spin frequency /(£) as a 
function of time (in d). The blue curve represents f(t) for (a) B = 2.28 X 1CT 8 , E = 9.28 X lCT 19 , p„ = 0.1, K = 1, Q n0 = 0,38 and 
fl = 0.68, and (b) B = 2.52 X 10" 8 , E = 1.16 X 10" 21 , p n = 0.1, K = 50, fl n0 = 0.65 and U = 0.84. The iMcCulloch et al.l j 1987ft 
data f bs(t) are graphed in red. The approximate solution (I22D is given by the dashed curve. The dotted curve is the steady-state spin 
frequency /(oo). 




p„(l+K) 

Figure 5. Superfluidity coefficients BE^ 1 ' 2 and p n (l + K) for all 
Vela glitches with N ^ 2. The bottom (blue) and top (red) lines 
correspond to I I26H as it stands, and with £i and ti swapped, re- 
spectively. The points correspond to by-eye fits to events with 
N > 3 (filled circles for Mt Pleasant, open boxes for Hartc- 
beesthock). 



longer present; there is just one combination of BE^ 1 ' 2 and 
p n (l + K) which leads to an acceptable fit. It is a testament 
to the power of the model that (i) it handles the triple-time- 
scale situation without introducing extra parameters, and 
(ii) the resulting fit constrains the superfluidity coefficients 
uniquely. 



3.4 Superfluidity coefficients 

The fitting procedure applied to the 1985 and 1988 test 
cases in £13.31 can be extended to all the events in Table [l] 
with JV ^ 2. In every instance, the two-component model 
matches the data at least as well as in Figures [2] and [4] The 
values of BE -1 ' 2 and pniX + K) from each event are plot- 
ted in Figure [S] The 16 pairs of diagonal lines come from 
plotting (|26|) for every glitch with N ^ 2, taking the time- 
scales t\ and ti both ways around; for glitches with N > 2, 
the two longest time-scales are assigned to £i and £2. The 
lower lines (blue) correspond to (|26[) . and the upper lines 
(red) correspond to swapping t\ and £2, as discussed in t]2.4l 



Measurements with different telescopes are plotted together. 
In addition, there are four events in Table [l] (1988, 1991, 
2000, 2004) with N > 3, for which the fitting procedure in 
SJ2J2 yields unique values of BE~ 1/2 and p n (l + K). The 
four independent measurements of these values are plotted 
as points in Figure [5] (open boxes for Mt Pleasant, filled 
circles for Hartebeesthoek). In all the fits, the very short 
(~ 1 min) time-scale £4 is excluded, as it is probably asso- 
ciated with diff e rent ( superconducting) physics; see [|6]and 
ISiderv fc Arparl (|2009D . 

Given that constitutive properties like viscosity and mu- 
tual friction should not vary significantly during the inter- 
val between glitches, we expect the values of BE^ 1 ' 2 and 
p n (l + K) inferred from all the observed Vela glitches to 
cluster. Figure 5 allows us to test this hypothesis. We find 
that the upper lines group around BE~ 1 ' 2 p~ 1 (l + A') -1 w 
87 ± 113 (mean ± standard deviation), while the lower lines 
group around BE~ xl2 p~ x (\ + K)' 1 « 0.268 ± 0.289. More- 
over, the four glitches with JV ^ 3, for which it is possible to 
determine BE -1 ' 2 and Pn(f + K) uniquely, group around 
BS _1/2 pn 1 (l + KY 1 ~ 0.80 ± 0.90. It is noteworthy that 
the unique fits favour the ordering £1 > £2; this is consistent 
with predictions from nuclear theory, as we shall explain in 

For a canonical light crust, with 0.3 ^ /On(l + K) ^ 3, 
the above results imply 26 ^ BE^ 1 ' 2 ^ 261 (upper lines) or 
0.08 < BE~ 1/2 < 0.8 (lower lines). In both cases, the mutual 
friction time-scale is within roughly one order of magnitude 
of the Ekman time-scale. We compare the above results with 
theoretical calculations of B and E in a 1 So neutron super- 
fluid and exotic quark matter in ijS] 



4 THE CRAB 

4.1 Data 

Since its discovery in 1968, the Crab pulsar has been moni- 
tored daily for glitch activity by a numbe r of groups (|Grothl 
ll975l ; lGullahorn et al.lll977l ; lLohsenlll98ll ). The most exten- 
sive and ongoing effort is overseen by the Jodrell Bank Ob- 
servatory, using a 12.5-m dish to measure all four Stokes 
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Table 2. Timing parameters for large glitches (Au ^ 1/xHz) in the Crab pulsar (i/ 9 o = 29.9 Hz). An asterisk in the CM column denotes continuous monit< 



Glitch Date 



1 
2 
3 
■1 
5 
6 
7 
8 
9 

10 
11 
12 
13 



MJD CM 



*3 



(1 



Av 3 



Av2 Av\ 

/xHz 



Ai/„ 



30-Scp-1969 

04-Feb-1975 

??-???-1981 

22-Aug-1986 

29-Aug-1989 

21-Nov-1992 

30-Oct-1995 

25-Jun-1996 

ll-Jan-1997 

10-Feb-1996 

30-Dcc-1996 

Ol-Oct-1999 

16-Jul-2000 



40494 
42448 
44900 
46664 
47767 
48947 
50021 
50260 
50459 
50489 
50813 
51452 
51741 



0.5 



18.7 
18 

9.3 
18 

2 
3.2 
10.3 
3.0 
2.2 
2.9 
3.4 
1.0 



97 
222 
123 

265 



-0.7 



-0.31 



Av Rcf. 



0.07 




0.05 


0.12 


1 


1.01 


-0.71 
-0.28 


1.02 


1.32 


1 
1 


0.12 


-0.11 


0.11 


0.12 


1 


2.28 


-2.11 


2.38 


1.85 


1 


0.26 




0.4 


0.3 


1 


0.064 




0.015 


0.08 


1 


0.66 




0.31 


0.66 


1 


0.2 




0.032 


0.23 


1 


-0.03 




0.05 


0.02 


1 


0.24 




0.017 


0.26 


1 


0.24 




0.04 


0.29 


1 


0.584 




0.143 


0.73 


2 



[1] (Wong et al. 2001) , [2] (Wang et al. 2001) 



parameters at 610 MHz for a maximum of 14 hr a day, 
occasionally supplemente d by readings at 1.4 GHz using 
the 76m Lovell telesco pe (|Lvne et al.ll 19951 : IShemar fc Lvnd 
ll996l ; lLyne et alj|200ch. A t otal of 26 glitches have been ob- 
served ( Melatos et al] 120081 ). Of these, 13 have been moni- 
tored sufficiently regularly to resolve the recovery stage. The 
timing parameters of the latter events are collected in Ta- 
ble [2] Although the 1969 and 1975 glitches were observed 
by several telescopes, Jodrell Bank was the only group to fit 
the data to the timing model given by (|17|l . Table[2]makes it 
clear that it is harder to resolve multiple, distinct, exponen- 
tial dec ays in the Crab th an in Vela. Only the 1989 glitch has 
N — 3 (|Lvne et al.lll992l ). and this event (the famous "slow 
glitch") is anomalous anyway, because the timing solution 
describes the spin-up event as well as the recovery stage, 
the only time the spin up has been resolved. As the spin up 
involves vortex physics outside the scope of the models in 
this paper, we exclude the 1989 glitch from our analysis. For 
all the other events, we have N ^ 2. 

Crab glitches differ from Vela glitches in a number of 
respects. The waiting-time probability distribution function 
is fitted accurately by an exponen tial in the Crab, i.e. it is 
consistent with a Poisson process (JMelatos et al] 120081 ). By 
contrast, Vela glitches quasi-periodically, wit h a small sub- 
set (~ 20%) of events spaced roughly evenly (|Melatos et al] 
I2008T ). The Crab glitches more frequently, at an average 
rate of 0.91^1 yr - 1 , compared to 0.43±S;i| yr _1 for Vela 
(jMelatos et al.ll2008f ). Its frequency jumps Au are measured 
to be ~ 1% of Vela's. During the recovery stage, the greatest 
distinction is the "overshoot" observed in the Crab, wherein 
the crust decelerates below its steady-state angular velocity 
before rising again asymptotically. This is characterized by 
the negative values of Av n in Table [2] which do not occur 
for Vela. The Crab has not attracted the same level of con- 
tinuous monitoring as Vela, so Table [2]is sparser than Table 
[T] More experiments are urgently required to investigate the 
form of the overshoot in more detail. 



4.2 Two-component superfluid 

We repeat the analysis in J)3.3l using the data in Table [2] To 
begin with, we apply the complete two-fluid model [equa- 
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Figure 7. Superfluidity coefficients BE -1 ' 2 and pn(l + K) in- 
ferred from the 1975 Crab glitch. The points are obtained by 
fitting the theoretical solution /(t) from (1 146 to the data by eye. 
The open boxes correspond to the fits in Figures [6f a) and (b). 
The filled circles indicate 12 additional fits. The two curves are 
plots of the relation (126 D obtained by fitting two exponentials to 
the data, with (ti,t 2 ) = (18, 97) (top) and (97, 18) (bottom) (tj 
and £2 in units of d). 



tion (|14|l ] to the 1975 glitch as a test case. In Figure [6] we 
present two distin ct fits to the timing solution measured by 
iLvne et al] (|1993h . Again, we find that the approximate so- 
lution (dashed curve) adheres more closely to the analytical 
solution (blue curve) in panel (b), where BE~ 1 ' 2 — 3.6 x 10 2 
is larger than in panel (a) (BE^ 1 ' 2 — 17). In both panels, 
we have p„ = 0.01 and K — 2500, demonstrating that for 
p n (l + K) = 25 there are two possible fits. These correspond 
to swapping the time scales £i and £2, as discussed in H3.3I 
Next, in Figure [7] we plot a range of fits obtained for 
the 1975 glitch on the BE~ 1/2 -p n (l + K) plane, c.f., Fig- 
ure^ for Vela. The parameters inferred by fitting the 1975 
glitch by eye are marked by solid circles, and the data from 
Figure [6] are marked by open boxes. Again, we find that 
(|26|l matches admirably the by-eye fits when there are two 
time-scales, even in the regime where the approximate solu- 
tion nominally breaks down. Indeed, the overall agreement 
is better than for Vela because the iterative numerical solu- 
tion to (|14[) is stable on the bottom branch too, cf. Figure 
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Figure 6. The 1975 Crab glitch, as fitted by a two-component spin-up model (see ij3.3H , showing the normalized spin frequency /(£) as 
a function of time (in d). The blue curve represents /(*) for (a) B = 7.3 X lCT 10 , E = 1.83 X 10~ 21 , p„ = 0.01, K = 2500, n n0 = 0,25 
and Q = 0.77, and (b) B = 3.3 X 10" 9 , E = 8.44 X 10" 23 , p n = 0.01, K = 2500, Q n0 = -2.47 and fl = 0.77. The lLvne et al.l jl993t l 
data f bs(t) are graphed in red. The approximate solution (I22D is given by the dashed curve. The dotted curve is the steady-state spin 
frequency /(oo). 



4.3 Overshoot 

Figure [6] makes it plain that the analytical model captures 
the "overshoot" in the Crab's recovery naturally. We can 
exploit the approximate solution (|22p to pin down the phys- 
ical conditions required for the overshoot to occur. As noted 
above, an overshoot requires Av n < for one value of n. In 
the context of equation ()22p . this translates into 



C 



fio — finO 



7BE- 1 / 2 /(20p n K) - 1 
or else 



<0, 



l-C-/(oo) = 



7BE-^ 2 {1 - Q )/(20p n K) + fi w0 - 1 
7BE~ 1 / 2 /(20p n K) - 1 



(27) 



< 0.(28) 



We can now distinguish four cases, depending on the sign of 
the denominator and whether (|27|l or (|28p is satisfied. If the 
denominator is negative, i.e. for (7BE~ 1 ' 2 )/(20p„K) < 1, 
and (|28[) is satisfied, we must have 1 < fio < finO- How- 
ever, this ordering produces a different type of overshoot, 
where the angular velocity of the crust initially increases 
above its steady-state value before decreasing asymptoti- 
cally. On the other hand, if (|27p is satisfied, the condition 
for an overshoot becomes fi„o < fio < 1- This is the case in 
Figure [Ha), where we have (7BE- 1/2 )/(20p„K) = 0.186, 
fino = 0.25 and fio = 0.77. Similarly, if the denominator 
is positive, i.e. for (7BE~ 1 ' 2 )/(20p n K) > 1, then we must 
have 1 < fio < fi„o if (|27l) is satisfied (which is not the 
overshoot we are looking for) and fi n o < fio < 1 if l|28p is 
satisfied. The latter is the situation in Figure [6jb), where 
we have (7BE- 1/2 )/(20p n K) = 5.39, fi„o = -2.47 and 
fio = 0.77. Therefore the condition for an overshoot of the 
form observed in the Crab is always 



fino < fio < 1 



(29) 



We can understand (|29l) physically in the context of 
Figure [6] In panel (a), the crust is light and it initially de- 
celerates rapidly below its steady-state angular velocity, in 
response to the viscous torque. Then, over a longer time- 
scale, mutual friction spins up the crust and viscous compo- 
nent into corotation with the inviscid component. In panel 
(b), the initial differential rotation between the viscous and 
inviscid components is rapidly removed by mutual friction. 




p„(l+K) 

Figure 8. Superfluidity coefficients BE~ X I 2 and pn(l + K) for 
all Crab glitches with N ^ 2. The bottom (blue) and top (red) 
lines correspond to (126 I t as it stands, and with t± and t,2 swapped 
respectively. 



While this is occurring, the even larger initial differential 
rotation between the viscous fluid and the crust drags the 
latter below its steady-state angular velocity. The viscous 
and invicid components achieve corotation first, followed by 
the crust, which is brought into corotation over the longer, 
viscous time-scale. Of these two possibilities, (a) appears to 
be the most natural; in the case of (b) a large and negative 
fi„o is required to generate an overshoot. 

The condition (|29p for an overshoot lends some insight 
into the glitch trigger. Firstly, glitches leading to an over- 
shoot can occur for all values of BE~ 1 ' 2 p„ 1 {K + 1) _1 in 
principle. Secondly, the internal fluid components must sat- 
isfy (|29p at t = 0+, i.e., immediately after the impulsive 
spin-up of the crust. The fact that overshoots are observed 
in the Crab suggests that (|29p is satisfied for the majority 
of its glitches, whereas ([29} is never satisfied in Vela. These 
differences must arise from the microphysics in the two ob- 
jects, such as the distribution of vortex pinning strengths 
which is plausibl y a function of age [e.g. if la ttice defects 
anneal over time (jMelatos fc Warszawskill2009l )]. 
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Figure 9. The 1975 Crab glitch, as fitted by a locked two- 
component spin-up model with B ~ 1 and p n = 0.01, showing 
the normalized spin frequency f(t) as a function of time (in d). 
The observational data / 6 s (i) are displayed as a heavy black 
curve. The lighter curves correspond to the theoretical model 
with K = 10 -2 (blue, top), K = 1 (green, centre) and K = 10 2 
(red, bottom). The dotted curve is the steady-state spin frequency 
/(oo). 



4.4 Superfluidity coefficients 

Finally, in Figure [8] we collect the inferred parameters for 
all Crab glitches with N ^ 2 on the BE~ 1/2 -p n (l + K) 
plane, c.f., Figure [5] for Vela. The four pairs of diagonal 
lines come from plotting (|26l) for every glitch with N = 2, 
taking the time-scales t\ and ti both ways around. The lower 
lines (blue) correspond to (|26|) as it stands, while the upper 
lines (red) correspond to swapping t\ and ti in ()26|) . We 
find BE- 1/2 p- 1 (l + K)- 1 « 38.5 ± 17.9 (mean ± standard 
deviation) for the upper lines, and BE~ 1 ' 2 p~ 1 (l + A') -1 fa 
0.270 ± 0.177 for the lower lines. For a canonical light crust, 
with 0.3 < p n (l+K) < 3, we therefore have 12 sC BE~ 1/2 sC 
116 (upper lines) or 0.08 ^ BE' 1 ' 2 ^ 0.81 (lower lines). 

Encouragingly, the Crab results are similar to Vela, 
lending general support to the model. The concordance 
is especially significant at t\ and £2 differ significantly 
in the two objects. In the upper fits, the mean value of 
B£ _1/2 p^ 1 (l + A") -1 for the Crab is about half that of 
Vela, while the mean in both pulsars is roughly the same for 
the lower fits in both objects. In general, typical values of 
B and E individually are slightly lower in the Crab than in 
Vela. For example, upon comparing Figure[6]with Figure[2l 
we find 10" 10 < B < 10~ 9 and 10~ 23 < E < 10~ 19 in the 
Crab, as against B fa 10~ 8 and 10" 23 < E < 10" 21 for Vela. 
In (Jll we discuss these numbers in the context of recent the- 
oretical and experimental constraints obtained from nuclear 
theory and heavy-ion collider experiments. 



4.5 Magnetically coupled crust and core 
superfluid 

The hydrodynamic model presented here can be used to rep- 
resent and test a variety of neutron star pictures. One im- 
portant possibility, discussed in ijl] £|2.1l and £|2-2t is that the 
viscous and inviscid fluid components lock together via mag- 
netic forces (e.g. if the proton-electron plasma is strongly 
coupled to the neutron superfluid via entrainment). This 
regime corresponds to B ~ 1 in our model (see £|2.2[) . Under 



such conditions, equation (|A1[) for the azimuthal velocity 
component of the interior flow reduces to 

V4, = -rcu + / dr'e" ; +( T " T ')/ (/) + rfi e"+ T (30) 

Jo 

with 



w+ 



pnJ(r) 



uj- = 



(31) 
(32) 



Importantly, because of (|32fl . there is now only a single 
time-scale in (|30p . which resembles the equation for a single- 
component fluid [see H3.2l and (|B1[) ], 

In Figure[9] we attempt to fit the 1975 Crab glitch using 
(fT3)l and IpO). As we can see, for 10~ 2 < A < 10 2 , the theory 
produces a nearly exponential decay on a single time-scale. 
Therefore, strong coupling (B ~ 1) de-activates the time- 
scale on which the viscous and inviscid components come 
into co-rotation, effectively reducing the problem to a single 
fluid. As in £13.21 this limiting case does not contain enough 
freedom to fit the observational data. Most importantly, we 
find that it never leads to an overshoot. Therefore, to re- 
produce the glitch recovery in the Crab, we must leave both 
time-scales in the problem, with their ratio to be determined 
by observations. 

Similar arguments can be made for the coupling of the 
viscous fluid (proton-electron plasma) to the crust. Theory 
predicts this coupling to be strong (~ 1 s) if the spin up is 
a result of induced tension in the magnetic field l ines, or 
~ 30 s if it is a result of classic Ekman pumping (|Eassonl 
1 19791 ). This scenario corresponds to if ~ 1 in our model. 
If the only remaining component is the neutron superfluid, 
it then responds on the mutual friction time-scale (Bfi) -1 . 
Again, only one remaining time-scale remains with which to 
fit the data. 

We emphasize that the failure of the above strongly 
coupled regime to fit the data does not mean that the neu- 
tron condensate and the proton-electron plasma are weakly 
coupled. The star might consist of multip le superfluid com- 
pone nts [e.g. multiple moments of inertia (jAlpar et alj|l993i . 
ll996T )] or a significant uncharged inviscid component (non- 
ideal excitations like in helium II, see ij2}. Such scenarios are 
contained in our model by taking K ~ 1. 



5 DISSIPATIVE PROCESSES IN BULK 
NUCLEAR MATTER 

The transport coefficients of bulk nuclear matter at ~ 
MeV energies have not yet been measured in terres- 
trial experiments. O f the 19 coefficients identified by 
lAndersson fc Comerl (|2006i ) in their flux-conservative for- 
malism, only the shear viscos ity 77 has been measured in rel- 
ativis tic heavy-ion colliders (|Adler et al J 12003 ; lAdare et al.l 
120071 ). These experiments (~ 500 nucleons at ~ 10 2 GeV) are 
conducted far from the physical conditions in a neutron star 
(~ 10 57 nucleons at ~ lMeV). Even so, they have thrown 
up the interesting finding that r\ approaches the quantum 
lower bound rj/s = h/inks (where s is the specific entropy) 
inferred from the duality between an ti-de-Sitter and confor- 
mal field theories JAdler et al.ll2003l ). 

In this section, we use the theory in ij2]and the fits to the 
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Vela and Crab data in Sj3]and SjHto constrain two transport 
coefficients: the shear viscosity 77, and the mutual friction co- 
efficient B. We also constrain the density ratio of the viscous 
and inviscid components, and hence, indirectly, the super- 
fluidity transition temperature. As the idealized model in ij2] 
does not incorporate stratification, the data are interpreted 
in terms of uniform rj and B, effectively representing mass- 
weighted averages of the depth-dependent r/(r) and B(r) in 
the real star. The reader is cautioned to bear this in mind 
when interpreting the findings. For example, if r\ and B are 
very different in the inner and outer core, the mass-weighted 
average may be a poor approximation to both regions. 




5.1 Neutron-rich matter 

In the outer core [1.6 x 10 14 < p/gcm~ 3 < 3.9 x 10 14 ], 
the neutrons condense into a So superfluid. In this phase, 
it is surmised that 77 arises from ele ctron-electron scat- 
terin g outside the superfluid vortices (|Cutler fe Lindbloml 
Il987l ), while B arises fr o m electrons sc attering off vortex 
cores (|Alpar et al.1 Il984l ; iMendelll Il99ll ). The bulk trans- 
port coefficients can be related to microscopic quanti- 
ties, like the charged-fluid-neutron-vortex relaxation time 
and Kelvin-wave oscillation frequency, in the low- frequency, 
long-waveleng t h lim i t by ge n eralis ing the method of 
lHall fe Vinenl (|l956l ). IMendelll (|l99ll ) calculated theoreti- 
cally that one has 



/j = 6.0 x lO^gcm'V 1 



10 14 gcm- 3 



B = 1.1 x 10 



B' = B 2 



7/6 



-2 (2/-!) x 
y 1 / 2 (l-x) 



uo 7 kJ 

1/6 



10 14 gcm- 3 



(33) 



(34) 



(35) 



where x is the proton fraction p n /p, and 0.3 ^ y — m* v jm v ^ 
0.7 is the normalized effective mass of the proton from Fermi 
liquid theory. For a cold equation of state below neutron 
drip, x is given by Eqn. (2.5.16) in IShapiro fe Teukolskvl 
1 19831 ). viz. 



with 



= ifi + 4 



4Q 4 (Q 2 - ml) 



-1 3/2 



+ 



(36) 



6.1 x 10 15 gcnr 



1/3 



(37) 



Q = m n — m p , and hence 2.6 x 10~ 3 ^ x/(l — x) ^ 0.125; 
the lower bound occurs at p s = 7.8 x 10 gem -3 . Clearly, 77, 
B, and B' are functions of depth through p, T, and hence x 
(and m*, weakly). 

In Figure [TCT1 we consolidate the ranges for BE" 1 ' 2 and 
p n (l + K) for Vela (Figure [5| and the Crab (Figure [8]) onto 
one diagram. The lightly shaded regions in Figure [lOl cover 
the range of Vela fits presented in Figure [S] In the lower re- 
gion, £1 and ti are given by (|24[) and (|25[) respectively; in the 
upper region, t\ and £2 are swapped. Similarly, the darkly 
shaded regions cover the range of Crab fits presented in Fig- 
ure [8] Overplotted are three theoretical curves for BE" 1 ' 2 
versus p n (l + K), constructed from equations (|33[) - (|37|) . To 



Figure 10. Ranges of the superfluidity coefficients BE ' and 
p n (l + K) predicted by nuclear theory [equations H33l l- H37l l1 and 
inferred from pulsar observations (Figures \E\ and |8j. The lightly 
and darkly shaded regions correspond to Vela and Crab data re- 
spectively. The theoretical curves are obtained by varying the 
density p for three different values of temperature, T = 10 7 ' 5 K 
(top, red), T = 10 70 K (centre, green) and T = 10 65 K (bottom, 
blue), with K = 50. 
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Figure 11. Theoretical mutual friction and viscous time-scales 
ti (blue, upper) and t^ < ti (red, lower), given by i l241 l and H25II 
respectively, as functions of pulsar spin-down age r c (in yr). In- 
put quantities E,p n and B are determined from nuclear theory 
by (1331) — <I37P >. assuming standard cooling to get a T-t c relation. 
Curves correspond to densities p = 10 14 (dashed curve), 10 15 
(solid curve), and 10 16 g cm -3 (dotted curve). The vertical bars 
denote the ranges of t\ and £2 measured observationally for the 
Crab (left) and Vela (right). 



plot these curves, we assume a canonical fluid-crust ratio, 
K — 50. (It is easy to slide the curves from left to right as 
K increases.) The three curves correspond to three different 
values of temperature, T = 10 75 K (top, red), T = 10 70 K 
(centre, purple) and T = 10 6,5 K (bottom, blue). The density 



10 12 gem" 



at 



10 17 gem 3 at the right-hand 



p s ~ p increases along each curve, from p 

the left-hand vertex, to p 

vertex. 

Figure [10] contains several important lessons. First and 
foremost, the three theoretical curves predicted by standard 
nuclear theory fall naturally on top of the parameter ranges 
inferred from the Crab and Vela data. This is strong circum- 
stantial evidence in favor of th e spin-down model in SjHand 
IVan Evsden fe Melatos! (|2010T ). It motivates further testing 
of the model as fresh data become available in the future. 
Second, the theoretical curves overlap with the shaded bands 
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in the regime BE~ 1 ^ 2 p~ 1 (l + A') -1 < 1, where the vis- 
cous time-scale £2 given by (|25[) is less than the mutual fric- 
tion time-scale £1 given by (|24p (|Mendelllll99ll : iReiseneggerl 
1 19931 ). Interestingly, the by-eye fits to the Vela glitches with 
N — 3 (points in Figure [5J also favour the £1 > £2 regime. 
Third, if we turn the argument around and treat 77 and B 
as known and given by (|33[) - (|37[) . then the intersection be- 
tween the theory and data constrains the fluid-crust ratio 
to 0.08 < Pn(l + K) < 4. The latter range is independently 
consistent with the standard nuclear equation of state. We 
discuss this issue further in £15-41 

The evolution of the recovery time-scales £1 and £2 with 
pulsar spin-down age r c can also be predicted from (1331) 
(|37[) . given a relation between T and r c . We assume standard 
cooling via the Urea and mo dified Urea p rocesses and a two- 
zone, heat blanket model (jPagel 1 199a ; iMelatos fc Peraltal 
120071 ; |Peraltall2007f ) . In Figure [Til we plot £1 and £2, defined 
by H24|) and (|25[). versus r c , using the theoretical formulas for 
viscosity and mutual friction in (f33)) - (f3"7)) and taking £1 > £2, 
in keeping with the conclusions of the previous paragraph. 
Results for three different values of p are presented; p = 
10 14 gcm -3 (dashed curve), p — 10 15 gcm~ 3 (solid curve) 
and p = 10 16 gcm~ 3 (dotted curve). The upper (blue) curves 
for £1 are flat (assuming m*/m p depends weakly on T). The 
lower (red) curves indicate that £2 decreases as T drops and 
77 decreases. Overplotted are the observed ranges of £1 and £2 
for Vela and the Crab, given in Tables [TJ and [2] respectively. 
We see that the theoretical curves are consistent with the ob- 
servations for 10 14 gcm~ 3 < p < 10 15 gcm~ 3 , nicely brack- 
eting the expected density of the outer core (jPeralta et al.l 
120051 ). The predicted decline of £2 for r c ^ 10 6 yr, where neu- 
trino cooling gives way to photon cooling, is an interesting 
test of the spin-down model. Glitches observed in pulsars 
older than ~ 10 6 yr are predicted to recover rapidly (< 1 d) 
immediately after the glitch (due to increased viscosity), fol- 
lowed by a slower recovery (> 10 d) mediated by mutual 
friction. 



5.2 Strange quark matter 



The analysis in £15.11 can be repeated for the numerous 
exotic fluid phases involving strange quark matter, which 
have been hypothesised to exist in neutron star interiors. 
Some of these phases are confined to the inner core, where 
the Ekman process may not penetrate due to stratif ication 
JAbnev fc Epsteinlll99d ; lvan Evsden fc Melatosll2008l '). Nev- 
ertheless, the analysis in this paper can distinguish between 
these phases in principle, by extracting values of 77 and B 
from glitch recovery data. It is therefore an aid in test- 
ing conjectures regardi ng the susceptibility of quark stars 
to r-m ode instabilities (IMadsen||2000l; iM annarcll i fc Manuel! 
120091 ') and starquake glitches JMannarelli et alj|2007l ')" 

In its simplest manifestation, deconfined strange quark 
matter is ungapped, i.e., diquark pairing does not occur. 
Under these circumstances, the shear viscosity arises pre- 
dominantly from quark-quark s cattering (modified by Lan- 
dau damping) and is given by drleiselberg fc Pcthick 1993; 
IJaikumar etaL 2008; Sa'd 2008; S hternin fc YakovlevlfaX)! ) 



V = 



\0.1J \n 



14/9 



T 



10 9 K 



-5/3 



-1 -1 
gem s 



(38) 



0.15 fm 3 is the nuclear saturation density (with p/m n ~ 
5 x no typically). For densities and temperatures pertinent 
to the core, e.g. 10 14 gem -3 < p < 10 15 gem" 3 and 10 8 ' 5 K < 



, e.g. 1U gem < p < 1U 1 

T < 10 95 K, equation (O predicts 4 x 10~ 15 < E < 6 x 

10- 13 . 

If diquark pairing does occur, the shear viscosity 
is suppressed exponentially by a fac tor exp( — A/3 fc_BT), 
where A is the pairing energy gap (Madscn 2000). For 
A > 1 MeV, the suppression is extreme. Several condensed 
phases have been identified, two of which we consider here: 
the color-flavor locked (CFL) and two-flavor superconduct- 
ing (2SC) phases. In the CFL phase (n < 10n , T < 
O.lMeV), u, d and s quarks pair up to form a conden- 
sate that is antisymmetric in it s color and flavor indices 
l| Afford et al.lll999l ; I Alfordl l200iT ) . Under these conditions, 
77 is dominated by phonon-phonon and kaon-kaon scatter- 
ing JMadsenll200ol; IJaikumar et al]|200a ; lAlford et ai]|2009l ; 
iMannarelli fc Manuell 120091 '); for T < O.lMeV, the elec- 
tron population is expo nentially suppressed and sc atters 
proportionately weakly (|Raiagopal fc Wilczekl 120011 ), For 
the phonon-phonon process (ignoring vortex scattering and 
"dr essed photons"), the shear viscosity is n ominally given 
by (|Manuel et al.ll2005l ; IJaikumar et al'j|2008l ) 



'-"**(&)' hkr 



-1 -1 

gem s 



(39) 



where p q « 300 MeV is the quark chemical potential. For 
temperatures and densities pertinent to the core, we obtain 
5 x 10~ n < E < 5 x 10~ 5 . Superficially, this phase appears 
to be ruled out by the consolidated data in Figure [TOl but 
it is vital to note that the phonon-phonon mean free path 
exceeds the stellar radius for T < 10 10 K, so the above value 
of 77 is only relevant for very hot cores, e.g. due to reheat- 
ing by r-modes. For the kaon-kaon process, assuming a toy 
interaction Lagrangian with a single kaon coupling constant 
(instead of the usual three), it is predicted that 77 depends 
microscopically on the Goldstone kaon phase speed and on 
whether the scattering occurs via a contact p rocess or the 
excha nge of a virtual particle. From Fig. 2 in lAlford et al.l 
(|2009j), we see that the kaon-kaon viscosity is between ~ 10 12 
and ~ 10 15 times lower than the phonon-phonon viscosity 
given by (|39[) . since the kaon chemical potential (« 4 MeV) 
is much less than p q . Importantly, however, for a typical 
core temperature T ~ 10 8 K, the phonon-phonon mean free 
path exceeds the stellar radius, so the kaons are the leading 
residual source of viscosity, with 10 -25 < E < 10~ 16 . The 
upper end of this range is consistent with the Crab and Vela 
dat a in Figure [ TOl 

IMannarelli et al.l (|2008bl ) showed that, for elastic 
phonon-vortex scattering (inelastic is suppressed) in the 
CFL phase, the mutual friction depends of the speed of 
sound, c B ~ 3 _1 ' 2 c, and the phonon density, with 



B = 



405c2 



c 2 



k B T 
"9 



(40) 



In (I38[) , a s denotes the QCD coupling constant and no 



For densities and temperatures in the core (see above), (|40l) 
evaluates to 5 x 10~ 21 < B < 5 x 10~ 16 . This is several orders 
of magnitude smaller than the typical B value extracted 
from the model in £13.41 and £14.41 

In the 2SC phase, the u and d quarks pair up to 
form a condensate and two quark colors acquire a gap 
(A ~ 50 MeV). However, one quark color remains ungapped, 
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Figure 12. The three curves are obtained by varying the density 
p for three different values of temperature, T = 10 7,5 K (top, 
red), T = 10 70 K (centre, green) and T = 10 6 ' 5 K (bottom, blue). 
Squares represent fits with K = 50 and while circles are for K = 1. 
The Vela fits arc shown in red, and the Crab is shown in blue. 



and this color inevitably dominates all dissipative processes 
in the system; the contribution to transport by the gapped 
colors is exponentially suppressed <x exp(— A/3/cbT), as in 
the CFL phase. (As no global symmetries are broken, the 
2SC phase is also not a superfluid.) Detailed enumeration 
of the interact i on channels available to the ungapped quark 
(|Madsenl I200Q ; Sa'd 2008) indicates that the shear viscos- 
ity of the 2SC phase is (5/9) ' 3 times the shear viscos- 
ity for standard, ungapped quark matter (with no diquark 
pairing), i.e., r\ is (5/9) 1 ' 3 times equation (|38l) . For densi- 
ties and temperatures in the core (see above), this implies 
4 x 10~ 15 <£<5x 10~ 13 . Our conclusions with respect to 
the glitch data are therefore the same, namely that the CFL 
phase (like standard, ungapped quark matter) is broadly 
consistent with the observations, especially at the lower end 
of the range. Ab initio theoretical calculations of B in the 
2SC phase have not yet been published, to the best of our 
knowledge. 



5.3 Depth averages 

An interesting property of the Crab and Vela fits is that they 
naturally yield values of the product BE -1 ' 2 and p n (1 + K) 
which agree nicely with nuclear theory, e.g. in Figures [lOl and 
1121 yet B and E disagree individually with nuclear predic- 
tions. One possibility, of course, is that the agreement in 
BE -1 ' 2 and p n (1 + K) is accidental; that is, the results are 
telling us that a self-consistent but purely hydrodynamic 
model is incapable of explaining the relaxation data, and 
new, non-hydrodynamic physics must be introduced. If true, 
this would fulfill one of the papers chief aims, as discussed 
in ijl] However, given that the discrepancies in B and E are 
not specific to this particular relaxation model but are the 
generic consequences of a hydrodynamic treatment, there is 
another, equally likely explanation: stratification and com- 
pressibility restrict Ekman p umping to a thin s u rface layer 
[not treated in this paper; cf . lAbnev fc Epstein! (| 19961 ) and 
Ivan Eysden fc Melatosl (J2008I VJ . so that the effective, body- 
averaged values of B and E used in the model are very dif- 
ferent to what one would expect for a uniform-density star. 
Figures [10] and [12] amply demonstrate this: the products 
BE -1 ' 2 and p„ (K + 1) depend weakly on density, varying 



~ 30-fold over the range 10 12 < p/gcnT 3 sC 10 17 (Fig- 
ure llOp . but individually they vary ~ 10°-fold over the same 
range of p. Furthermore, to bring the B and E fits into closer 
agreement with nuclear theory, one would need to shift the 
neutron-rich-matter curves in Figure[l2]to the left and down, 
just as one expects if the Ekman layer is restricted. 

Figure \T%\ displays the fitted values of B and E from 
jj3] and if4] for all Crab (blue points) and Vela (red points) 
glitches, assuming K — 1 (circles) and K = 50 (squares). 
Overall, both B and E are smaller than the theoretical pre- 
dictions in (|34l) <I36[) for neu tron-rich ma t ter. T he coupling 
time t„ ~ 1 s determined by lAlpar et al.l (| 19841 ) for typical 
values of density, proton density fraction and effective pro- 
ton mass translates to B ~ 3 x 10~ 4 and B ~ 7 x 10" 



using Mendell's relation [see Eq (21) in iMendelll (|l99ll )] 
for the Crab and Vela respect ively. S i milarl y, the respec- 
tive Ekman times estimated by lEassonl (|1979l ) are translate 
to E ~ 10 -6 and 10 -5 , marked as blue and red crosses 
in Figure 1121 Compositional variations such as the exis- 
tence of strange quark matter in the interior also affect the 
depth-averaged effective parameters, as discussed in H2.ll 
The ranges of B and E for the CFL phase from H5.2I are 
sketched on figure [12] as a dashed rectangle. They are un- 
certain, of course, but depth-averaging over strange quark 
matter dramatically reduces the effective B and E. Figure 
1121 illustrates that the self-consistent hydrodynamic model 
potentially possesses considerable discriminatory power be- 
tween nuclear models which are widely separated on the 
B-E plane, once it is refined further. 



5.4 Crust fraction 

The theoretical curves in Figures [10] and [11] are drawn for 
the canonical value K — 50 of the crust-fluid fraction. 
How does this compare with terrestrial experiments that 
probe the equation of state of neutron-rich nuclear matter? 
iLattimer fc Prakasbl (J2007T ) showed that K depends sensi- 
tively on the symmetry energy at subsaturation densities 
through the transition pressure p+ and density p+ at the 
solidification point. Quantitatively, one has 



1 
A" 



28^> +J R 3 (l-1.67£-0.6£ 2 



3M»c 2 



£ 



1 + 



2p+ (1 + 5C - 14g 2 
p + m b c 2 t; 2 



with £ = GM t /R*c 2 , where mj, M*, and R* denote the 
mean baryon m ass, stellar mass , and stellar radius respec- 
tively. Recently, IXu et al.l (|2009n analysed isospin diffusion 
data from heavy-ion collisions and placed bounds on the den- 
sity dependence of the symmetry energy. They concluded 
that p+ and p+ lie in the ranges 0.01 MeVfm -3 ^ p+ ^ 
0.26MeVfm~ 3 and 0.040 fm" 3 < p+/m n < 0.065 Me Vfm" 3 , 
consistent with t he isotopic de pendence of giant monopole 



resonances in S n (ILi et al. 12007 ) and the neutron-skin thick- 
ness of 208 Pb (|Li fc Chenl 120051 ). The above ranges imply 
86 $J K ^ 1779. Again, this prediction is consistent with 
the results in Figure 1101 allowing for moderate shifts of the 
theoretical curves to the right. It is also co nsistent with the 
bound K < 70 inferred from Vela glitches (|Link et al.lll999T ) 
invoking angular momentum conservation. 



(41) 



Pulsar glitch recovery and the superfluidity coefficients of bulk nuclear matter 15 



6 CONCLUSIONS 

A two-fluid hydrodynamic model which predicts the spin 
evolution of the crust during the recovery stage of a pulsar 
glitch is presented. The model has six free parameters (p n , 
K, B, E, fi,jo and Qq) and solves for the back reaction of 
the viscous torque on the crust in a self-consistent manner. 
This analysis repr esents the ne x t step in the class of mod- 
els introduced by ISiderv et al.l IJ2009T ) , by replacing body- 
averaged internal components with an exact global Ekman 
flow pattern. 

The paper presents a straightforward recipe for extract- 
ing two important transport coefficients in bulk nuclear mat- 
ter, namely the shear viscosity and mutual friction parame- 
ter, from radio timing observations of glitch recovery. Given 
a post-glitch timing solution of the multi-exponential form 
(|17p . or equivalently (|19|) . one fits the data to the theoreti- 
cal expression (|22[) and reads off four of the model param- 
eters (e.g. B, E, $l„o and Qq) in terms of the other two 
(e.g., K, p n ) from equations igDJ, IJ2JJI and ^-^. The 
fits (e.g. in Figures [2] [4] and [6| are excellent and include 
a natural explanation of the "overshoot" phenomenon seen 
in the Crab. For glitches with three resolvable time-scales, 
the model parameters are constrained uniquely, and sat- 
isfy (7BE~ 1/2 )/(20p n K) < 1. This implies that the viscous 
time-scale, appearing in (|25|) . is shorter than the mutual fric- 
tion time-scale, appearing in (|24|) . When the fits derived for 
all glitches with N ^ 2 in Vela (Figure[S]) and the Crab (Fig- 
ure |HJ are plotted together on the BE~ 1 ' 2 -p n (K + 1) plane, 
the regions covered by the two objects overlap. Moreover, 
the region of overlap coincides with curves of BE -1 ' 2 versus 
p n (K + 1) predicted from nuclear theory, both for a stan- 
dard So neutron superfluid at the temperatures and densi- 
ties expected in the outer core and also for certain phases 
of strange quark matter (e.g. stable ungapped phase, 2SC 
phase, or CFL phase with kaon viscosity) under inner core 
conditions. The crust-fluid ratio inferred from the same fits 
is consistent with rece nt measurement s of isospin diffusion 
in heavy- ion collisions (|Xu et al.ll2009f ). 

The self-consistent spin-down model in this paper helps 
solve two of the puzzles raised in iJTJ namely that glitch 
recovery cannot be parameterized by a simple exponential 
decay and that it is not always monotonic. ft does not solve 
the third puzzle - that the recovery time-scale varies from 
glitch to glitch in the same object - but it does shed new 
light upon it by sharpening the debate. Of the six parame- 
ters, the transport coefficients B and E do not change sig- 
nificantly during the interval between glitches. The inertia- 
related factors p n and K may conceivably change with each 
event, if different zones within the star (e.g. hosting dif- 
ferent a verage pinning strength s ) participate in different 
glitch es IjAlpar et al.l Il993l . 1 19961 ; ISedrakian fc Hairapetianl 
|2002j). However, the global nature of the Ekman process ar- 
gues against this; even if successive glitches are triggered in 
distinct pinning zones, the resulting global flow embraces 
the entire neutron condensate and charged fluid (and exerts 
a torque on the whole crust), so that the effective values of 
p n and K are always the same. This leaves one possibility: 
fino and fio vary from glitch to glitch. Such an outcome is 
natural if different pinning zones trigger different glitches; 
the shear preceding a glitch is not always the same (e.g. we 
do not observe a reservoir effect in the data). 



The model presented in this paper possesses several se- 
rious weaknesses. First, it fails to explain the fourth time- 
scale £4 ~ min measured by iDodson et al.l (|2002i . 120071 ) 
in two Vela glitches. This phenomenon is probably a re- 
sult of physics not included in this analysis, e.g., the re- 
laxation of magnetic flux tubes (|van Hoven fe Levinl 120081 ; 
ISiderv fc Alparl [2009J) , or a solid core in a superconduct- 
ing p hase. Second, the two-fluid theory neglects entrain- 
ment (|Andersson fc Comerll2006MSiderv et alfeoooh . to keep 
the difficult spin-down problem tractable analytically, even 
though entrainment is known to be important. Our conclu- 
sions regarding r\ and B are expected to change by factors 
of a few due to this effect. Third, the analysis neglects bulk 
viscosity, because it enters at order O(E) in the equations of 
motion and therefore does not modify the Ekman process. 
In practice, however, the bulk viscosity exceeds the shear 
viscosity by several orders of magnitude in several phases of 
neutron-rich and exotic quark matter as when electroweak 
processe s like d + s o u + s work to restore chemical equi- 
librium jMadsenll200(il ; lManuel et al.ll2005l ; Ijaikumar et all 
l200Sl ; ISa'dll2008l ). Lastly, although the values of BE~ 1/2 ex- 
tracted when fitting the model to data are consistent with 
what nuclear theories predict, the individual values of B and 
E are not. The most natural explanation is that B and E 
depend sensitively on p whereas the combination BE^ 1 ' 2 
does not, and that stratification and compressibility effects 
limit the Ekman flow to a narrow surface layer in a way that 
is not modeled properly in this paper. This important issue 
deserves further study. 

More radio timing data on glitch recoveries is urgently 
required to fill out Figure [TO] and test the theory more 
comprehensively. Multibeam monitoring experiments with 
phased radio arrays are under consideration and offer excit- 
ing prospects in this regard. 
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APPENDIX A: GENERAL SOLUTION TO THE 
TWO-FLUID EKMAN PROBLEM 

The general solution to fTJ-Q in a n arbitrary container for 
any ch oice of B, B' , E, p n is given bv lVan Eysden fc Melatosl 
(|20ld ) 



v^(r,r) 



rui+(r)oj-(r) 



dr'f(r') 



(Al) 



/3[w+(r) -w-(r,j J0 

x |[o; + (r) + /3]e"+ M (— ') 

rQ n0 uj + (r)ij-(r) r M+ ( r ) T _ w_(r)T-| 

uj-{r)e +y J — uj+(r)e y J 



P[u+(r)-u)-(r)] 



w+(r) -w_(r) 
where we have defined 
2B£T 1/2 



B> 



(A2) 



which characterises the relative strength of the superfluid 
mutual friction force to the viscous forces, and the exponen- 
tial time-scales 



W±(r) = -- 



P- 



m 

h(r) 



*ti 



P + 



h(r) 



PpnJ(r) 
h(r) 



.(A3) 



which are a combination of the classical Ekman time-scales 
and the mutual friction time-scale. In (IA3I). we have 



I(r) = 



X.(r)[X 2 _(r) + X 2 + (r)] 



(A4) 



x{ 



J(r) = 



1 - H\r) 

i + m{r) 

H\r) 



[A 2 + (r)-A 2 _(r)] 2 +4A^(r)A 2 _(r)} 1 / 2 

(A5) 



2A_(r)[l + # 4 (r)] 

x {[A 2 _(r) + A 2 + (r)] [l - H 4 {r)] + 4X1 (r) H 4 (r)} 



and 






A±(r) = 


H(r)\ 


(B'-2) 2 + B 2 


_{p n B' - 2) 2 + {p n B) 2 






p s B(l + H 4 (r)) 



1/2 



1/2 



(A6) 



H 2 (r)[(p n B'~2) 2 + (p n B) 2 ] j 

The factors H(r) and h(f) are related to the shape of the 
boundary. In a general container which is symmetric about 
the z — plane, at the upper boundary z = h(r) while at 
the lower boundary z — —h(r). Then H(r) is defined as 



(r) = {l+[ft'(r)] 2 } 



H(r) = 



1/4 



(A7) 



For a sphere, they are given by (|Van Eysden &, Melatosl 
l2010h 



h(r) = H{r)~ 2 = (l-r 2 ) 



1/2 



(A8) 



The viscous torque on the crust is computed from the 
right hand side of 



df(r) 
dr 



15K 



a 2, , sdv<t,(r) 
drr h(r) 



<:h 



(A9) 



Equations (|A1[) and (|A9[) . combine to give an integral equa- 
tion for the spin evolution of the crust, viz. 

/(r) = -p n K f dr' [g A (r - r') + g B (r - r')] f{r') (A10) 



+ p n K [g A (r)Q n0 + g B {r)n ] + 1 : 



where 



At \ 

9 (t) 



15 

2 



1 r 3 J(r) r e "+M T - e "-M-l 
dr — 



o w-(r)-w+(r) 

g B {r) = p I dr'aV). 



(All) 



(A12) 



In the limit B,B' < 1 (but not necessarily BE 1/2 < 
1), we find that (|A2|) - (|A6|) simplify to 

P = BE' 1 ' 2 (A13) 

I{r) = J{r) = X- ± 1 {r) = H{r) (A14) 

which give the results (fTTj) - (fT6)) . 



APPENDIX B: GENERAL SOLUTION OF THE 
ONE-FLUID EKMAN PROBLEM 

When there is no superfluid component, we have B — B' = 
0, p n = 1, n„o = fio- Therefore I(r) = J(r) = Aj x (r) = 
H(r), P = 0, uj+ = 0, and oj_ = —H(r)/h(r). Equation 
(IA1I) then reduces to 
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r jm / ^H^-m^ 1 ). 



Mr,r) = -rff dr7(r')e 'H l ; +rS! e »M(B1) 

^y) Jo 

which is the c l assica l solution obtained by 
I Greenspan fc Howard! (| 19631 ). In a sphere, the bound- 
ary is defined by (JMJ). Equations ()A10|| - (|A12|) then 
become 



f(r) = -K / dr'.g A (r-r')/(r') (B2) 

Jo 

+ K 9 a (t)Q + 1, 



g A (r) = — - j drr A h(r) 



.mil, 

e Hr) 



(B3) 



g B (r) = 0. (B4) 

The only variables in this limit are K, fio, and E. Equations 
((20)) and ([21]) simplify to 



(B5) 



Ai/p _ 1 + KQ.0 
T^J ~ 1 + K ' 

*y N , Ais n /t n 20 i io 

providing two conditions in terms of the observed quantities 
on the left-hand sides. We are left with one model variable 
(say, K) which is then determined by fitting the shape of 
/(r) to fobs{r) by eye. 



APPENDIX C: PARAMETERS EXTRACTED 
FROM BY-EYE FITS 

Table [Cl quotes the model parameters p n , K, B, E, Qq an d 
Qno for all glitches fitted by-eye. 

This paper has been typeset from a TpX/ I^Tf?X file prepared 
by the author. 
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Table CI. By-eye model fits to glitch data 



Object Glitch Ref. 



A 



B x 10" 



E x 10" 



Op n n o 



Crab 



1975 



0.1 


80 


8.5 


1.30 


0.77 


0.127 


0.1 


500 


7.0 


0.049 


0.77 


0.28 


0.1 


750 


6.5 


0.024 


0.77 


0.31 


0.1 


1000 


6.36 


0.014 


0.77 


0.32 


0.1 


1500 


6.36 


0.0062 


0.77 


0.32 


0.1 


2500 


6.36 


0.0023 


0.77 


0.33 


0.01 


2500 


7.3 


0.18 


0.77 


0.25 


0.1 


80 


32 


0.083 


0.77 


-2.46 


0.1 


500 


33.5 


0.0020 


0.77 


-2.55 


0.1 


750 


34 


0.00087 


0.77 


-2.59 


0.1 


1000 


31 


0.00049 


0.77 


-2.59 


0.1 


1500 


31 


0.00022 


0.77 


-2.59 


0.1 


2500 


31 


0.00008 


0.77 


-2.59 


0.01 


2500 


33 


0.0084 


0.77 


-2.47 


0.1 


1 


228 


92.8 


0.68 


0.38 


0.1 


50 


252 


0.116 


0.84 


0.65 


0.001 


500 


253 


14.7 


0.84 


0.69 


0.001 


2000 


253 


0.78 


0.84 


0.67 


0.001 


5000 


253 


0.13 


0.84 


0.67 


0.01 


2000 


253 


0.0078 


0.84 


0.67 


0.01 


5000 


253 


0.0013 


0.84 


0.67 



Vela 



1985 



Vela 



1988 



0.01 



53 



50 



3.05 xlO 5 



0.97 0.97 



[1] (Lvne et al. 1993) . [2] (McCulloch et al. 1987) . [3] (Flanagan 1990) 



